Hive Parameter Analysis

library(tidyverse)
library(ggplot2)
library(readr)
library(viridis)

Loading required package: viridisLite

library(ggfortify)


round_any <- function(x, accuracy, f=round){f(x/ accuracy) * accuracy}

bee_param_df <- list.files(path="hive_data/heat_queen/", full.names = TRUE) %>% 
                lapply(read_csv, show_col_types = FALSE) %>% 
                bind_rows

New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:New names:

write.csv(bee_param_df,"bee_data_isaac_heat_queen.csv",row.names = FALSE)

What these all mean: n = queen cells per hour rb = brood radius rn = necter radius w = total daily honey pph = pollen ratio ph = honey consumption ratio pp = pollen consumption ratio k = consumption probability value

broodMetric = average number of brood surrounding other brood pollenRing = average min distance between honey and brood

brood_bee_param <- read.csv("bee_data_isaac_heat_brood.csv") %>% mutate(type = "Brood") #%>% mutate(wha = 0, qha = 0)
worker_bee_param <- read.csv("bee_data_isaac_heat_worker.csv") %>% mutate(type = "Worker") #%>% mutate(bhd = 0, qha = 0)
queen_bee_param <- read.csv("bee_data_isaac_heat_queen.csv") %>% mutate(type = "Queen") #%>% mutate(bhd = 0, wha = 0)
none_bee_param <- read.csv("bee_data_isaac_heat_none.csv") %>% mutate(type = "None") #%>% mutate(bhd = 0,wha = 0, qha = 0)

bee_heat_param_df <- rbind(brood_bee_param,worker_bee_param) %>% rbind(.,queen_bee_param) %>% rbind(.,none_bee_param)

ggplot(bee_heat_param_df, aes(x = type, y = pBroodHeat))+
  geom_boxplot()+
  theme_classic()

heat_glm <- glm(pBroodHeat ~ 0 + type, data = bee_heat_param_df)
heat_glm

Call:  glm(formula = pBroodHeat ~ 0 + type, data = bee_heat_param_df)

Coefficients:
 typeBrood    typeNone   typeQueen  typeWorker  
   0.07445     0.22406     0.06848     0.16027  

Degrees of Freedom: 944 Total (i.e. Null);  940 Residual
Null Deviance:      37.67 
Residual Deviance: 14.31    AIC: -1266

Looking at correlation between heat parameter and brood in the place

rand_bee_param <- read.csv("bee_data_isaac_rand.csv") %>% mutate(type = "Rand")

Warning: cannot open file 'bee_data_isaac_rand.csv': No such file or directoryError in file(file, "rt") : cannot open the connection

bee_param_df_perc <- bee_param_df %>% select(c(type,pBrood,pHoney,pPollen,pEmpty)) %>%
                                      pivot_longer(c(pBrood,pHoney,pPollen,pEmpty))

ggplot(bee_param_df, aes(x = broodMetric, y = pollenRing, color = type))+
  geom_point(alpha = 0.25)+
  ylim(0,16)+
  xlim(0,6)+
  theme_classic()

ggplot(bee_param_df_perc, aes(x = type, y = value, fill = name))+
  geom_boxplot()+
  theme_classic()

ggplot(m1_bee_param, aes(x = pPollen, y = pph))+
  geom_point()+
  theme_classic()

broodMetricGLM_m0 <- glm(broodMetric ~ 1, data = m2_bee_param)
broodMetricGLM_m1 <- glm(broodMetric ~ n + rb + rn + w + pph + ph + pp + k, data = m2_bee_param)
broodMetricGLM_step <- step(broodMetricGLM_m0, direction = "both", scope = formula(broodMetricGLM_m1), trace = 0)

summary(broodMetricGLM_step)

Call:
glm(formula = broodMetric ~ n + rb + k, data = m2_bee_param)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-0.85986  -0.11668   0.02039   0.13922   0.48804  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.2085892  0.0629850  66.819   <2e-16 ***
n            0.0075472  0.0005096  14.810   <2e-16 ***
rb          -0.0256635  0.0109140  -2.351   0.0191 *  
k            0.0040269  0.0020931   1.924   0.0550 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.0380154)

    Null deviance: 26.924  on 481  degrees of freedom
Residual deviance: 18.171  on 478  degrees of freedom
AIC: -202.19

Number of Fisher Scoring iterations: 2

pollenRingGLM_m0 <- glm(pollenRing ~ 1, data = m2_bee_param)
pollenRingGLM_m1 <- glm(pollenRing ~ n + rb + rn + w + pph + ph + pp + k, data = m2_bee_param)
pollenRingGLM_step <- step(pollenRingGLM_m0, direction = "both", scope = formula(pollenRingGLM_m1), trace = 0)

summary(pollenRingGLM_step)

Call:
glm(formula = pollenRing ~ k + rn + pph + rb + n + w + ph, data = m2_bee_param)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  
-2.23594  -0.24853   0.02459   0.27886   1.69306  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  9.915e+00  4.021e-01  24.656  < 2e-16 ***
k            9.135e-02  5.177e-03  17.646  < 2e-16 ***
rn          -3.120e-01  1.983e-02 -15.738  < 2e-16 ***
pph         -8.050e-01  9.764e-02  -8.245 1.64e-15 ***
rb          -2.297e-01  2.706e-02  -8.488 2.72e-16 ***
n           -8.865e-03  1.265e-03  -7.006 8.46e-12 ***
w           -1.369e-04  2.719e-05  -5.036 6.76e-07 ***
ph          -9.297e-01  3.702e-01  -2.512   0.0123 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for gaussian family taken to be 0.2321617)

    Null deviance: 292.08  on 481  degrees of freedom
Residual deviance: 110.04  on 474  degrees of freedom
AIC: 673.91

Number of Fisher Scoring iterations: 2

library(viridis)

Loading required package: viridisLite

ggplot(queen_pos_df, aes(x = QueenX, y = QueenY))+
  geom_path()+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df, aes(x = QueenX, y = QueenY, color = Dist2Heat))+
  geom_point()+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df, aes(x = QueenX, y = QueenY, color = cos(AngleOppHeat)))+
  geom_point()+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df, aes(x = QueenX, y = QueenY, color = sin(AngleOppHeat)))+
  geom_point()+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df_heatmap, aes(x = QueenXRound, y = QueenYRound, color = cos(meanAngleMoved)))+
  geom_point(size = 5)+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df_heatmap, aes(x = QueenXRound, y = QueenYRound, color = sin(meanAngleMoved)))+
  geom_point(size = 5)+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

ggplot(queen_pos_df_heatmap, aes(x = QueenXRound, y = QueenYRound, color = count))+
  geom_point(size = 5)+
  geom_point(aes(x=11,y=48.5),colour="red") +
  scale_color_viridis() +
  theme_classic()

hive_pca <- prcomp(bee_heat_param_df %>% select(-c(...1, trial_n, type, days, n,
                                                   pBrood, pPollen, pHoney, pEmpty,
                                                   pBroodHeat, pPollenHeat, pHoneyHeat, pEmptyHeat,
                                                   broodMetric,pollenRing)),
                   center = TRUE,scale. = TRUE)
summary(hive_pca)

Importance of components:
                          PC1    PC2    PC3    PC4    PC5     PC6     PC7     PC8     PC9    PC10
Standard deviation     1.0604 1.0538 1.0513 1.0167 1.0039 0.99266 0.98240 0.96346 0.95014 0.91471
Proportion of Variance 0.1124 0.1111 0.1105 0.1034 0.1008 0.09854 0.09651 0.09283 0.09028 0.08367
Cumulative Proportion  0.1124 0.2235 0.3340 0.4374 0.5382 0.63672 0.73323 0.82605 0.91633 1.00000

autoplot(hive_pca, colour = "pHoney", loadings = TRUE, loadings.label = TRUE,
         data = bee_heat_param_df)+
  scale_color_viridis()+
  theme_classic()

ggplot(bee_heat_param_df, aes(pPollenHeat, pBroodHeat))+
  geom_point()+
  theme_classic()

ggplot(bee_heat_param_df, aes(wha, pPollen))+
  geom_point()+
  theme_classic()